# Make graphs for figure 1

# Load Models and Store AUCs
countries <- c('colo', 'indo')
dvs <- c("any", "high", "spike")
titles <- c(any = "Any violent event",
            high = "Five or more violent events",
            spike = "1SD Increase in events")

for (country in countries) {
  for (v in dvs) {
    source(paste("config_code/",country,"_data_setup.R",sep = ""))
    source(paste("config_code/",country,"_predictor_vars.R",sep = ""))
    models <- c(names(figures_1_2)[-1], names(figure_A5))
    T <- end.year - start.year + 1
    
    # Get population benchmark dataset
    pop_auc_array <- array(NA,
                           dim = T)
    filename <- paste(modeldir,"/",
                      country,
                      "_ebma_",
                      v,
                      "_population.RData",
                      sep = "")
    load(filename)
    for (i in 1:T) {
      predictions <- as.vector(ebma.results[[i]]$fit.oos)
      realizations <- as.vector(ebma.results[[i]]$actual.oos)
      pop_auc_array[i] <- roc(response = realizations,
                              predictor = predictions)[['auc']]
    }
    
    auc_array <- array(NA,
                       dim = c(length(models),
                               T),
                       dimnames = list(models,
                                       NULL))
    for (x in models) {
      filename <- paste(modeldir,"/",
                        country,
                        "_ebma_",
                        v,
                        "_",x,".RData",
                        sep = "")
      load(filename)
      for (i in 1:T) {
        predictions <- as.vector(ebma.results[[i]]$fit.oos)
        realizations <- as.vector(ebma.results[[i]]$actual.oos)
        auc_array[x, i] <- (roc(response = realizations,
                                predictor = predictions)[['auc']]
                            - pop_auc_array[i])
      }
    }
    mean_auc <- apply(auc_array, 1, mean)
    
    # Make Graph
    fixed <- names(figures_1_2_fixed)[-1][!(names(figures_1_2_fixed)[-1]
                                            %in% shock_breakdown)]
    slow <- names(figures_1_2_slow)[!(names(figures_1_2_slow)
                                      %in% shock_breakdown)]
    annual <- names(figures_1_2_annual)[!(names(figures_1_2_annual)
                                        %in% shock_breakdown)]
    var.order1 <- order(mean_auc[fixed],
                        decreasing = FALSE)
    var.order2 <- order(mean_auc[slow],
                        decreasing = FALSE)
    var.order3 <- order(mean_auc[annual],
                        decreasing = FALSE)
    n_breakdown = length(shock_breakdown)
    var.order4 <- n_breakdown:1
    
    ranks1 = rep(NA,length(var.order1))
    ranks1[var.order1] <- ((1:length(var.order1))
                           +length(var.order2)
                           +length(var.order3) + n_breakdown
                           +ifelse(length(var.order2)>0,
                                  3,2))
    ranks2 = rep(NA,length(var.order2))
    ranks2[var.order2] <- (1:length(var.order2)
                           +length(var.order3)
                           +n_breakdown
                           +2)
    ranks3 = rep(NA,length(var.order3))
    ranks3[var.order3] <- 1:length(var.order3)+n_breakdown + 1
    ranks4 <- n_breakdown:1
    ranks <- c(ranks1,ranks2,ranks3, ranks4)
    varsets.total <- c(fixed,
                       slow,
                       annual,
                       shock_breakdown)
    pdf(paste("graphs/",
              country,
              "_AUCs_",
              v,
              "_ebma_popplus_ShockBreakdown.pdf",
              sep=""),
        width = 7, 
        height = 7)
    par(mar=c(5, 12, 2, 1) + 0.1)
    plot(mean_auc[varsets.total],
         ranks,
         xlim = c(-.1,.1),
         ylim = c(0,max(ranks)+1),
         ylab = "",
         xlab = "Performance improvement over population (AUC)",
         yaxt='n',
         main = "")
    for (i in 1:length(ranks)) {
      lines(x = c(-1,1),
            y = c(ranks[i],ranks[i]),
            col = "gray95")
    }
    points(as.vector(t(auc_array[varsets.total,])),
           rep(ranks,each = T),
           pch = rep(c(17,
                       rep(16,T-2),
                       15),
                     length(varsets.total)),
           cex = rep(c(1,
                       rep(.5,T-2),
                       1),
                     length(varsets.total)),
           col = rep(c("gray",
                       rep("black",T-2),
                       "gray"),
                     length(varsets.total)))
    points(mean_auc[varsets.total],
           ranks)
    axis(2,
         at = ranks,
         labels = labels_A5[varsets.total],
         las = 2)
    axis(2,
         at = c(length(ranks)+ifelse(length(var.order2)>0,
                                     4,3),
                n_breakdown +length(ranks3)+2,
                n_breakdown + 1,
                n_breakdown + length(ranks2)+length(ranks3)+3
         ),
         labels = c("Time-Invariant",
                    "Time varying (Annual)",
                    "Shock Breakdown",
                    ifelse(length(ranks2)>0,"Slow-Moving","")),
         las = 2,
         font=4 )
    lines(x = c(0,0),
          y = c(-1,length(ranks)+ifelse(length(ranks2)>0,7,5)),
          lty = 3)
    dev.off()
  }
}
